home *** CD-ROM | disk | FTP | other *** search
/ Eagles Nest BBS 8 / Eagles_Nest_Mac_Collection_Disc_8.TOAST / Developer Tools⁄Additions / MacScheme20 / Mathlib / roots.scm < prev    next >
Encoding:
Text File  |  1989-04-27  |  2.6 KB  |  68 lines  |  [TEXT/????]
  1. ;;; The following procedure attempts to improve on the precision of the
  2. ;;;  roots returned by POLY->ROOTS, and will adequately find multiple
  3. ;;;  roots as well. (POLY->ROOTS appears to split these multiple roots
  4. ;;;  by perturbations on the order or (sqrt *machine-epsilon*); in some
  5. ;;;  contexts this splitting may be a desirable feature.) It does this
  6. ;;;  by computing the first two derivatives yp and ypp of the given
  7. ;;;  polynomial y, and using a single newton-iteration on the function
  8. ;;;  u = y/yp (which is why ypp is required). We use u instead of y
  9. ;;;  because it handles the problem of multiple roots: such a root of y
  10. ;;;  becomes a single root of u. This is essential because only a single
  11. ;;;  newton-iteration is to be used.
  12.  
  13. (define (new-improved-poly->roots y . initial-guess)
  14.   (let* ((fuzzy-roots (apply poly->roots (cons y initial-guess)))
  15.          (yp (deriv-poly y))
  16.          (ypp (deriv-poly yp)))
  17.     (define (newton x)
  18.       (let* ((yx (poly-value y x))
  19.              (ypx (poly-value yp x))
  20.              (yppx (poly-value ypp x))
  21.              (ux (/ yx ypx))
  22.              (upx (- 1 (* yx yppx (/ 1 ypx ypx)))))
  23.         (- x (/ ux upx))))
  24.     (map-n newton *fuzzy-repeat-count* fuzzy-roots)))
  25.  
  26. (define (new-improved-poly->roots y . initial-guess)
  27.   (let* ((fuzzy-roots (apply poly->roots (cons y initial-guess)))
  28.          (yp (deriv-poly y))
  29.          (ypp (deriv-poly yp)))
  30.     (define (newton x)
  31.       (let* ((yx (poly-value y x))
  32.              (ypx (poly-value yp x))
  33.              (yppx (poly-value ypp x))
  34.              (ux (/ yx ypx))
  35.              (upx (- 1 (* yx yppx (/ 1 ypx ypx)))))
  36.         (- x (/ ux upx))))
  37.     (map-n newton *fuzzy-repeat-count* fuzzy-roots)))
  38.  
  39. (define (nu-poly->roots y . initial-guess)
  40.   (let* ((fuzzy-roots (apply poly->roots (cons y initial-guess)))
  41.          (yp (deriv-poly y))
  42.          (ypp (deriv-poly yp))
  43.          (num (mul-polys y yp))
  44.          (den (sub-polys (mul-polys yp yp) (mul-polys y ypp)))
  45.          (qr (div-polys num den))
  46.          (q (car qr))
  47.          (r (cadr qr)))
  48.     (define (newton x)
  49.       (- x (+ (poly-value q x)
  50.               (/ (poly-value r x) (poly-value den x)))))
  51.     (map-n newton *fuzzy-repeat-count* fuzzy-roots)))
  52.  
  53. (define *fuzzy-repeat-count* 1)
  54.  
  55. (define (improved-poly->roots y . initial-guess)
  56.   (let* ((fuzzy-roots (apply poly->roots (cons y initial-guess)))
  57.          (yp (deriv-poly y)))
  58.     (define (newton x)
  59.       (let ((yx (poly-value y x))
  60.             (ypx (poly-value yp x)))
  61.         (- x (/ yx ypx))))
  62.     (map-n newton *fuzzy-repeat-count* fuzzy-roots)))
  63.  
  64. (define (map-n f n list)
  65.   (if (zero? n)
  66.       list
  67.       (map-n f (- n 1) (map f list))))
  68.